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ABSTRACT 



Context. Energetic gamma rays (GeV to TeV photon energy) have been detected toward several supernova remnants 
(SNR) that are associated with molecular clouds. If the gamma rays are produced mainly by hadronic processes rather 
than leptonic processes like bremsstrahlung, then the flux of energetic cosmic ray nuclei (> 1 GeV) required to produce 
the gamma rays can be inferred at the site where the particles are accelerated in SNR shocks. It is of great interest to 
understand the acceleration of the cosmic rays of lower energy (< 1 GeV) that accompany the energetic component. 
These particles of lower energy are most effective in ionizing interstellar gas, which leaves an observable imprint on the 
interstellar ion chemistry. A correlation of energetic gamma radiation with enhanced interstellar ionization can thus 
be used to support the hadronic origin of the gamma rays and to constrain the acceleration of ionizing cosmic rays in 
SNR. 

Aims. We propose a method to test the hadronic origin of GeV gamma rays from SNRs associated with a molecular 
cloud. 

Methods. We use observational gamma ray data for each SNR known to be associated with a molecular cloud, modeling 
the observations to obtain the underlying proton spectrum under the assumption that the gamma rays are produced 
by pion decay. Assuming that the acceleration mechanism does not only produce high energy protons, but also low 
energy protons, this proton spectrum at the source is then used to calculate the ionization rate of the molecular cloud. 
Ionized molecular hydrogen triggers a chemical network forming molecular ions. The relaxation of these ions results in 
characteristic line emission, which can be predicted. 

Results. We show that the predicted ionization rate for at least two objects is more than an order of magnitude above 
Galactic average for molecular clouds, hinting at an enhanced formation rate of molecular ions. There will be interesting 
opportunities to measure crucial molecular ions in the infrared and submillimeter-wave parts of the spectrum. 



Key words. Astroparticle physics - Radiation mechanisms: non-thermal 
supernova remnants - Gamma rays: ISM 
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1. Introduction 

The origin of cosmic rays (CRs) is an open question in as- 
trophysics. The cosmic ray spectrum below the knee, at en- 
ergy/nucleon E < 10^^ eV is believed to be associated with 
cosmic ray acceleration in supernova remnants (SNRs) (Bell 
T978bl IBlandford fc Ostriker| [19781 [Blandford fc Eichier 



19871. However, there is no conclusive proof for this un- 
til now. There is an excess in GcV-TeV gamma rays ob- 
served from SNRs associa ted with mol ecul ar clouds, see e.g. 



Abdo et al. ( 2009 J , Abdo et al. ( 2010c ) and Aharonian et al 
( |2008p . These signals might be caused by bremsstrahlung 
or inverse Compton scattering of electrons in a leptonic sce- 
nario, or by the decay of neutral pions formed by proton- 
proton scattering in a hadronic scenario. So far, it is not 
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known which of these processes is dominant. Investigating 
which is the dominant process is important to understand 
the origin of cosmic rays. In the hadronic scenario, high en- 
ergy protons are accelerated in the SNR shocks and then 
escape to interact with ambient protons, in particular in 
molecular clouds in the direct vicinity of the SNR. It is also 
likely that low energy protons {E < 1 GeV) are accelerated 
in the SNR, but these protons fall below the threshold for 
pion formation, so that no conclusions concerning the low 
energy CR spectrum can be drawn from gamma ray ob- 
servations. However, low energy protons are very efBcient 
in ionizing molecular gas. Therefore, ionization signatures 
provide information about the density of low energy cos- 
mic rays. The main product of the ionization of molecu- 
lar hydrogen, H^, initiates a chain of chemical reactions 
that yield additional ions like H^, 0H+, H2O+, H3O+ and 
HeH+ ( ,Black,2007j McCarthy et al.|2006[ ) . These molecules 
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are most likely formed in rotatioiially and vibrationally ex- 
cited states, the corresponding wavelength for relaxation 
in the UV or IR, respectively. If the abundances of these 
molecules are sufficiently large, the UV or IR signals might 
be detectable and offer conclusions concerning the source of 
cosmic rays. A correlation study of molecular clouds bright 
in GeV gamma rays and showing ionization features might 
be useful to find the dominant process in forming GeV 
gamma rays. In this paper, the ionization rate of molec- 
ular hydrogen is calculated for each SNR known to i nteract 



with a mol ecular cloud. In contrast to former work (Becker 
et al.|20lT I , here the spectral shape of the primary particle 



spectral energy distribution (SED) below kinetic energies 
of ~ 1 GeV is altered in order to take the unknown spec- 
tral behavior into consideration, rather than altering the 
minimum energy of particles contributing to the ionization 
process. The paper is structured as follows: In section [2] 
the competing processes active while particles are being 
accelerated and their influence on the spectral shape of the 
CRs is discussed, in section |3] the spectrum of the primary 
protons is calculated by considering loss processes and the 
acceleration mechanism, in section [4] the ionization rate for 
each SNR associated with a molecular cloud is calculated, 
in section [5] the uncertainties are discussed, in section |6] the 
ionization signatures to be expected are shown, in section 
[7] first observational evidence for an enhanced ionization 
rate in correlation with GeV gamma ray emission is sum- 
marized and in section [8j a summary of the paper as well 
as an outlook to future work is given. 



2. Acceleration and diffusion 

Since the primary particle spectra at energies below 
1 GeV at the source are not known, especial ly in the 



context of cosmic ray-induced ioniz ation (e.g. Nath & 
|Bier man"ir (1 994[ ) , Indriolo et al. ( 2009 ) ) , the competing pro- 
cesses affecting these particles are discussed and compared 
in this section. In particular, it is of importance at what 
energy the ionization timescale is shorter than the acceler- 
ation timescale and vice versa to ensure that acceleration 
is una ffected. T he acceleration timescale is given in |Jokipil] 
( 1987[ ), [Biermann et al.| ( |2009| ) as 



8k 



(1) 



where k is the diffusion coefficient and Vsh is the shock 
velocity. The diffusion coefficient may well differ inside the 
cloud and outside the cloud. Outside the cloud, the diffusion 
coefficient has to be low enough to allow for efficient acceler- 
ation, while inside the cloud the diffusion coefficient has to 
be sufficiently large for the particles to penetrate the cloud 
within the age of the SNR. For a typical age of T = 10"^ yr 



and a penetration depth of i? = 30 pc ( Becker et al.|20lT ) 
the diffusion coefficient for a particle of p = 

would have to be k = ~ lA x 10^* 



2 T 

Introducing a momentum dependence in the diffusion co- 
efficient, K = kq ( 1 GeV c-i ) ' ^^'^ applying this value 
for the cloud's interior, the acceleration timescale can be 
written as 



= 4.4 X 10 
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500 km s" 



P 



1 GeV c' 



(2) 
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Fig. 1. Ionization and acceleration timescale for protons, 
n = 100 cm^'^ and Veh = 500 km/s. 



where c is the speed of light and p is the momentum of the 
particle. This includes a scaling to a typical value of the 
shock velocity for m iddle-aged remnants, V^sh = 500 km 
(jAbdo et al.| j2010c |). 

The momentum loss rate dp/dt at a given momentum p 
by ion ization or excitation is given by |Lerche fc Schlickeiser 
(fT982| as 



dp 
di 



-5 X 10-19 



11.3 + 2 -In 



V cm 3 y V mc I 



\mc) 



(3) 



eV cm 



where q is the charge of the particle, n is the number density 
of the interacting region, m is the mass of the particle and 
c is the speed of light. A more general ex pression can be 
found in Mannheim & Schlickeiser (1994). Neglecting the 
logarithmic dependence in the square bracket, since p is 
of the order of ~ 1 GeV 0"^, the ionization timescale for 
protons is calculated as 



•p = 5.6 X 10^^ 



V y 

mc) 



(4) 

These two timescales are compared in Fig. [TJ where g = 1, 
n — 100 cm"'^ and T4h = 500 km/s were used as a typical set 
of parameters. The momentum dependence of the diffusion 
coefficient is shown fo r 8 = 1/3 as well as for 8 = 0.6, as 
discussed iu jBlasi fc Amato (2012] and references therein. 
There it is reported that a value of o = 1/3 would favor sec- 
ond order Fermi acceleration and at the same time explain 
the observed ratios of B/G and fit the anisotropy of cosmic 
rays observed at Earth better than 8 — 0.6. However, this 
would require an injection spectrum of N{E) oc E~^ '^, 
which is challenging for non-linear diffusive shock accel- 
eration (NLDSA) in SNRs. On the other hand, 5 = 0.6 
would favor NLDSA and match the detected spectra of cos- 
mic rays including nuclei heavier than helium, but result 
in an anisotropy larger than observed. The shock velocity 
Vsh = 500 km/s i s typical for older SNRs, as e.g. W51C 



( |Abdo et al. 

can reach va 



20091. For younger SNRs, the shock velocity 
ues of Fsh ~ 10,000 km/s. 



In the environment described by the chosen parameters, 
at a particle momentum of p > 0.8 GeV c~^, almost in- 
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dependent of the actual momentum dependence of the dif- 
fusion coefficient, the acceleration timescale is shorter than 
the ionization timescale, indicating that ionization losses 
do occur, but do not suppress the acceleration process ef- 
fectively above this momentum. Therefore, the ionization 
losses do affect the spectral index of the primary proton 
SED significantly, but only at momenta p < 1 GeV c~^. 

Furthermore, adiabatic deceleration might in principle 
occur and alter the spectral shape of the prima ry proton 
SED, especially fo r momenta p> I GeV (see Lerche & 
Schlickeiser ( 1982 )). The momentum loss for primary parti- 
cles with momenta p > 1 GeV by adiabatic deceleration 
dominates the ionization losses, while below this energy the 
ioni zation losses dominate l osses by adiabatic deceleration 
(see Lerche & Schlickeiser (1982)). However, the primary 
proton SED at these energies is obtained from modeling 
the Ge V gamma ray emission via Tr'^-decay using the K amae 
model (Kamae et al. 2006 Karlsson & Kamae 2008). The 
lowest observable photon energy E^y « inD~Mev corre- 
sponds to a primary proton energy of Ep = 1 GeV. Thus, 
above 1 GeV the spectrum of primary protons needed is di- 
rectly known and would already include modulation effects 
on the SED caused by adiabatic deceleration. Yet, there is 
no direct information about the low energy cosmic rays, so 
the estimate of the low energy cosmic ray spectrum has to 
be made very carefully. 



3. Primary proton SED 

To calculate ionization by cosmic ray protons in a molec- 
ular cloud, the cosmic ray proton spectrum at the cloud 
is required. There is no observational data of the particle 
spectrum at the source: only the spectrum after propaga- 
tion to Earth is known. Due to larger uncertainties in the 
description of the propagation process, especially due to the 
lack of knowledge of the exact magnetic field configuration, 
the spectrum at the source cannot be described easily from 
the observed data. If the magneti c field is known, it can b e 
taken into consideration following 'Padovani & Galli (2011 ). 
Gamma ray emission from hadronic interactions is, on the 
other hand, very well suited to derive the primary parti- 
cle spectrum above ~ 1 GeV, since the gamma spectrum 
follows the primary spectrum. This gamma ray emission 
was detected e.g. by the Fermi-LAT instrument. Assuming 
that the detected gamma radiation is mainly caused by the 
decay of neutral pions from inelastic proton-proton inter- 
actions at the cloud, the primary proton SED can be found 
modeling the gamma ray detections. The spectral shape of 
this SED is gained modeling the gamma ray emission from 
neutral pion decay, induced by inelastic proton-proton in- 
teractions, and fitting the modeled gamma ray spectrum to 
the observational data. Loss processes for the primary par- 
ticles are summarized in subsection 
cesses are discussed in subsection 13.2 



decay, on the other hand, provides the primary proton spec- 
trum at the location of the formation of the neutral pions, 
the molecular cloud, since they decay in less than 10 
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3.1[ acceleration pro- 
and finally the calcu- 



( [Particle Data Group|2010| ). The gamma rays emitted from 
the decay of these pions do hardly suffer energy losses. 
Therefore, the primary proton spectrum at the location of 
the cloud is obtained. Low energy protons are very likely 
accelerated in the same place as the high energy protons, so 
in this cloud there is not only the formation of pions, but 
also ionization to be expected. Furthermore, deceleration 
of high energy proto ns additionally increas es the number of 



low energy protons (Padovani et al. 2009[ ). Since the spec- 
trum obtained from modeling holds for the location of the 
molecular cloud, no additional momentum losses have to be 
considered. However, the primary proton SED can only be 
considered as known from the modeling down to energies 
of roughly 1 GeV, as mentioned above. Below this energy, 
there is no observational information about the primary 
proton spectrum available. Ionization, on the other hand, 
will only be effective below 1 GeV, with the cross section for 
the direct ionization of molecular hydrogen by an incom- 
ing proton peaking at about 10 keV and rapidly dec lining 
with increasing proton energy (Padovani et al. (2009j) and 
references therein). This makes estimates for the primary 
proton spectrum below ~ 1 GeV rather uncertain and a 
crucial part of the calculation of the ionization rate. It is 
well possible and most probable that the power law behav- 
ior of the primary proton SED is not to be extrapolated 
to lower energies, because the acceleration mechanisms in 
these energy domains may differ. In order to account for loss 
processes which are effective for energies below ~ 1 GeV 
as well as the unknown acceleration mechanism at these 
energies, here the primary proton spectrum derived from 
modeling is attenuated to lower energies by the choice of a 
broken power law with positive spectral index 1 < a < 2, 
E^"", which is compatible with predictions of Skilling fc 



Strong ( 1976 ) , for three diff erent lower break energies Ei\^ 
It should be mentioned that Zirakashvili & Ptuskin ( 2008 ) 



lation of the normalization of the primary proton SED is 
described in detail in subsection 13.31 



and Ellison & Bykov (20II ) predict a concave spectrum for 
the primary particles escaping SNRs and entering nearby 
molecular clouds. The models used there apply for young 
SNRs with an age of 10^ yrs or younger, while here all 
examined SNRs are at least middle-aged, abo ut 10* yrs or 
older. For old SNRs, |Bozhokin fc Bykov] ( |1994[ ) modeled the 
penetration of a broad cosmic ray proton spectrum into a 
molecular cloud, assuming a diffusion coefficient with mo- 
mentum dependence k oc p^-^^ and k oc p^-^ . Their results 
motivate further examin ation of SNR s inter acting with a 
molecular cloud, as e.g. Bykov et al. (2000) modeled the 
spectrum of cosmic ray electrons interacting with a molec- 
ular cloud for the case of IC443, where they also derived 
a profile of the ionization rate due to pri mary elect r ons in 
the shocked part of the cloud. Recently, Yan et al. (2011) 
investigated the acceleration of protons in SNRs and the 
generation of gamma rays in nearby molecular clouds, tak- 
ing into account the streaming instability and background 
turbulence. 



3.1. Loss processes 



3.2. Acceleration mechanism 



The primary protons accelerated by the SNR can suffer 
momentum losses on their way to the molecular cloud. Yet, 
it should be stressed here that the spectral shape gained 
from modeling the gamma ray emission via neutral pion 



Fermi acceleration (Fermi 1949) is the very first approach 
to stochastically accelerate charged particles at shock fronts 
and the model has been further deve loped to the theory 
of diffusive shock acceleration (DSA) ( Krymskii||1977 Bell 
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1978b|a[ pandford fc OstrikCTl[l978l |Schlickeiser|[l989a|b ). 

A variety of concrete models concerning the accelera- 
tion mechanism at work in SNRs have been estabhshed 
the 



pa s t (e.g 



. Scott fc C hevalier (1975), Blandford 
|& Eichler' (I l987|"', l lilasil pOOSft, |Zirakashvih fc Aharonian| 
pOTO), Eich ler fc Pohi|(|201ip,|Drury|P011| ) aiIJre?e rences 



therein) . Recently, Uchiyama et ai. "( |2010p described an al- 
ternate model where the focus is not on escaped CR par- 
ticles, but strong adiabatic compression behind the SNR 
shock wave in the cloud leads to diffusive shock acceleration 
of pre-existing cosmic rays. With their model, they are ca- 
pable of modeling both the observed synchrotron radiation 
and gamma ray emission from certain SNRs associated with 
a molecular cloud, in particular W51C, W44 and IC443. 
However, there is no direct evidence which process actually 
is responsible for the acceleration. Furthermore, propaga- 
tion effects and possibly reacceleration are important. The 
unknown nature of the actual acceleration mechanism leads 
to a huge variety in possible CR spectra be l ow ^ 1 GeV, 



which is shown in figure 1 of Indriolo et al. (2009), where 



several propagated primary spectra are summarized. But 
since the modeling of the primary particle spectrum fitting 
the observed gamma ray emission via neutral pion decay 
does offer the spectral shape of the primaries at the source 
above 1 GeV, the primary spectra used in this work are in- 
dependent of the actual acceleration mechanism. The only 
assumption made is that the major contribution to the GeV 
gamma ray emission from the sources is caused by neutral 
pion decay, which shall be tested this way. 



Ey + dE-y) and a unit volume from the decay of neutral 
pions is reported as: 



dE^ dV dt 



«H / ainciiEp)jp{Ep)F^ 



E, 



dEp 
E„ 



where nn is the density of the ambient medium and 
CineiC^'p) is the cross section of inelastic proton-proton in- 
teraction^ The function Fj{x, Ep) describes the number 
of photons in the energy interval {x^x + dx) per colli- 
sion and is a dimensionless probability density distribution 
function. 

The result of the formula mentioned above is the num- 
ber of gamma rays formed from neutral pion decay per 
unit time, unit energy and unit volume at the location of 
the hadronic interactions, in units of GeV""'^ s~^ cm~'^. To 
transform this quantity into the quantity detected, the re- 
sult first has to be multiplied by the volume of the interac- 
tion region in order to obtain the total number of gamma 
rays formed from neutral pion decay per unit time and unit 
energy: 



dA^^ 
dE^ dt 



cloud 



(yinc\{Ep)jp{Ep)F^ 



E-y 



Er 



Ep J Ep 



dK, 



(9) 



3.3. Calculation of the proton 5ED normalization 

In this subsection a detailed description of the calculation 
of the normalization of the proton SED using the observed 
gamma spectrum is given. Known is the flux of gamma rays 
detected at Earth, 



^ dE^ dt d^Ea 



rth 



(5) 



in units of GeV~^ s^^ cm^^. What can be calculated from 
observations is the normalization of the flux of high energy 
protons interacting with the cloud and thus forming neutral 
pions causing the gamma radiation, which is continued to 
lower energies in order to calculate the ionization rate of the 
cloud next to the SNR. This flux is needed not at Earth but 
at the SNR, 



dA„ 



dEp dt dAs 



ap*p(^p) 



in units of GeV s cm . For instance, for W51C, 

-1.5 / rn \ -1.4 



$p(^p/ 



E^ 



1 GeV 



1 



15 GeV 



(6) 



(7) 



Here, $p(i?p) is a dimensionless spectral function and Op 
is the normalization factor in units of GeV~^ s~^ cm~^. 
The latter enters the calculation of the ionization rate. In 
order to obtain this value, the proton flux at the source is 
calculated from the observed gamma ray spectrum. 

The calculation of the formation of gamma rays via 
neutral pion decay is described in detail in e.g. |Kelner| 
fc Aharonian (2008). In this paper, the equation for the 
formation rate of gamma rays in the energy interval {E^^ 



Assuming that these gamma rays are emitted isotropically, 
the fraction that is detected at Earth per unit area can be 
gained from this by dividing by {'^''^d'^a.j:th.-ao\ircc) account 
for the solid angle: 



— ^'){E~^)Vc\oud 



(And 



Earth— source/ 



(10) 



Using eqs. (|8| and (10), the photon spectrum at 
Earth can be written as 



where 



<JineliEp)<S>p{Ep)F^ 



Op^HVcloud 



Er 



dEp 
E„ 



^" "Earth-source 



= const. 



(11) 



(12) 



is the normalization constant of the gamma spectrum. This 
factor is determined by high-energy gamma observations of 
the sources. 

For the calculation of ionization rates, the proton SED 
normalization, Op has to be calculated. This is done using 
conservation of energy: 

<PiEp) Ep dEp 



dNp 

dEp dt dylgource 



Ep dEp — ftp 



where v{Ep) = I 1 — ( 1 



Wp 

Kloud 
1/2 



(13) 

c is the velocity of 



the particle depending on its kinetic energy Ep and Wp is 



the CR fiux is used. 



Kelner & Aharonian ( 2008 \ use the CR density, while here 
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the total proton energy budget of protons with a minimum 
energy of E'min- Solving this for Up gives: 



HEp) 



(14) 



primary proton ionization is significant at the considered 
energies. The other ionization processes by protons have 
cross sections which are significantly lower than the one 
for direct ionization by primary protons, as can be seen in 
figure 1 of Padovani et al. (2009). 



where ap is in units of GeV~^ s~^ cm~^. So the entire 
expression for the normalization of the gamma spectrum 
from neutral pion decay, a^, can be written as 



^''^ "'Earth-source 



$(gp) Ep dEp 



= const. 



(15) 



The gamma normalization is therefore independent of the discussed in IPadovani et al 



Ionization by primary electrons is neglected here for 
two reasons: First, the ionization cross sections for these 
processes are lower than the corresponding ones for pro- 
tons. Second, while protons hardly lose energy on their way 
from the SNR to the molecular cloud, electrons can suf- 
fer energy losses. However, since the primary CR spectra 
near the peak of the corresponding ionization cross section 
(~ 10^ eV) is unknown, it is possible that primary elec- 
trons do con tribute to the ionizatio n rate significantly, as 

(2009). Because the focus of 



cloud volume, while the proton SED normalization, ap, de- ^^.^ ^^^^ on the contribut iSE^CR protons, here it is 



pends on this volume 

The model ing of the gamma rays is done using the 
Kamae model ( Kamae et al. 2006 Karlsson & Kamac 2008|) 



in order to obtain the value of a^. Here, a factor of 1.85 Ts 
multiplied to the resulting spectrum in order to take he- 
lium a nd heavier nuclei into account, as suggested by |Mori 
(2009). Here, the model is modified by these two factors 
and will be referred to as the modified Kamae model. 

For a given spectral shape of the proton SED, $p(£'p), 
the gamma ray spectrum normalization, a^, can be found 
from modeling. For a given distance of the object from 
Earth, dEarth-sourcej and lower integration threshold, i^niin, 
this leads to a certain value for the product Wp ■ nn for 
each object. Because there are no precise estimates of the 
average hydrogen densities of the objects, nn, here a value 
of TiH = 100 cm~^ is assumed. The result for Wp simply 
scales inversely with njj, if should turn out to be differ- 
ent. With a value for Wp, the proton SED normalization Up 
is calculated for each object. Since ap is already implicitly 
including the solid angle interval, the multiplication by 47r 
in equation (17) is not performed. 



4. Calculation of the ionization rate 

In general, ionization by particles is mainly caused by two 
different kinds of particles, namely electrons and protons. 
One has to consider direct ionization by electrons as well as 
protons on the one hand and ionization by electron capture 
of protons on the other hand. The full expres sion for the 
ionization rate following Padovani et al. ( 2009 ) is 



= 47r 



k 



jk{Ek)[l + MEk)]ar{Ek)dEk 



-4tt 



jp{Ep)a;-'^-{Ep)dEp , 



where -Emin is the minimum energy of particles consid- 
ered to contribute to the ionization process, jk is the num- 
ber of CR particles of species fc (A: = e or p, respectively) 
per unit time, area, solid angle and energy interval, (7)^°" is 
the ionization cross section for particles of species fc, CTp'^' 
is the electron capture cross section for protons and 4>k is a 
number taking into account that ionization may not only be 
due to ionization by a primary particle fc, but also by the 
electrons set free during this ionization, called secondary 
ionization. A closer look at the orders of magnitude for the 
different summands shows that only the contribution from 



assumed that the contribution of primary electrons to the 
total ionization rate is dominated by the contribution of 
primary protons and secondary electrons. In fact, ioniza- 
tion by secondary electrons is an important aspect, but it 
only increases the ionization rate by less than a factor of 
2. The resulting ionization rates derived here are rather 
lower limits, due to neglecting the contribution of primary 
electrons. 

Since only primary particles with a minimum kinetic 
energy of 10^ eV can penetrate the cloud (as will be dis- 
cussed below), this contribution seems negligible. On the 
other hand, protons penetrating the cloud will lose energy 
gradually, so in principle the effect of secondary ioniza- 
tion needs to b e calcu lated differentially, as was done in 
[Padovani et al. (2009). This will be done in future work 
and thus the ionization rates calculated here are lower lim- 
its. 

Neglecting the aforementioned ionization processes, the 
calculation reduces to: 



An 



jp{Ep)al°^iEp)dEp 



(17) 



For the examined SNRs the lower limit of integration, E'rnin, 
depends on the hydrogen density. This is due to the fact 
that the CR particles have to penetrate the molecular cloud 
before interacting with the gas, and lower energy particles 
are decelerated on shorter length scales than higher ener- 
getic ones and therefore do not contribute noticeably to the 
ionization. However, deceleration of high energy particles 
populates the low energy part of the spectrum. The stop- 
ping length of the particles depends on the hydrogen den- 



sity, an d thus does the lower integration limit. |Indriolo et al. 



(2009) suggest a lower integration thre shold of 2 MeV for 
diffuse clouds as originally suggested by |Spitzer fc Tomasko| 
(1968), riH < 10^ cm~^, and a lower integration threshold 



(16) of 10 MeV for dense clouds, tt-h > lO'^ cm ^. But since this 



effect is not well understood in quantitative terms, 100 MeV 
for dense clouds or 100 keV for diffuse clouds in the direct 
vicinity of the SNR might be considered as well. The direct 
io nization cross sec tion for primary protons used is reported 
by Padovani et al. (2009) in eqs. (5) and (6). 



Using the proton spectra obtained from modeling the 
gamma ray detections, the ionizati on r ate by primary pro- 
tons is calculated from equation (17 1. These spectra are 



varied in the spectral index a below the lower break en- 
ergy Eib. Additionally, the lower break energy Ei\, is var- 
ied. The minimum energy of protons penetrating the cloud 
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and thus contributing to the ionization process consid- 
ered is -Ernin = 10 McV as a conservative approximation. 
The Galactic average ionization rate of molecular hydro- 
gen in molecular clouds is C^^i aver ~ ^ xlO^^^ s^^, as 
reported by Neufeld et al. (2010). It should be noted that 
this average value is subject to variations between approxi- 
mately 10.5 xlO~^^ s~^ at maximum and 0.5 x 10~^^ s~^ or 



less, possibly connected to propagation effects (Indriolo & 
McCallJ|2012 ). The proton spectra of the sources are either 
of the form 



JpiEp) = 



or of the form 



[Ep < Sib) 

{Ep > Eih), 
(18) 



JpiEp) = 



«p(t)"^exp(-^)(|^) 

«p(S)''-p(-iat^) 



where Op is the normalization factor, 



En 



(Ep < Sib) 

{Ep > Sib), 
(19) 

= 1 GeV or 

1 TeV (see Table [l]), Sbr is the location of the spectral 
break, s = a; is the lower spectral index, As-|-s = is the 
higher spectral index, Scutoff is the higher cutoff energy, Sib 
is the location of the lower spectral break and a = 2.0, 1.5 or 
1.0 the spectral index below the lower break Sib- It should 
be mentioned that the spectral break in the primary proton 
spectrum is partly due to the fact that it is given in terms 
of the kinetic energy Ep. Diffusive shock acceleration pro- 
duces a power law in momentum, so expressing this spectral 
behavior in terms of the kinetic energy Sp, the spectrum 
deviates from a power law in kinetic energy near the rest 
energy of the particle, ~ 1 GeV for protons. Because the 
particles are accelerated at the supernova shock front which 
in some cases penetrates the molecular cloud, both acceler- 
ation and ionization can occur in the same place and at the 
same time. To account for the unknown acceleration mech- 
anism below S ~ 1 GeV, the lower spectral break Sib is 
introduced. Below this break energy, the particle spectrum 
is assumed to decrease rapidly toward lower particle ener- 
gies as to give a conservative lower limit on the ionization 
rate. For Sib < 1 GeV, this does not change the resulting 
gamma ray spectrum from neutral pion decay. 

Here, jp(Sp) = dA'^p/(dSp dt dAsomcc) is the number of 
CR protons per unit time, area and energy at the source. 
The spectral shape of the proton spectrum <I>p(Sp) for each 
object is found modeling the gamma ray detections by 
Fermi- L AT (see refer e nces [Abdo et al.| (|2009|), |Abd o et al. 
(|2010c), Abdo et all (|2010al, Abdo et al.| m 



^ .20lOc|[ ), [Abcfo 

[eTal., (2010b| ), |Castro & Slane (2010))" assuming hadronic 
interactions to form neutral pions, which decay via gamma- 
gamma coincidences to be the cause of the gamma ray emis- 
sion. The spectral information about all sources is given in 
Table El 
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To calculate the primary proton SED, observational 
data about the distance d of the object from Earth and 
the a ppr oxi mate volume of the cloud is required (see equa- 
and 



due to the maximum value of the ionization cross section 
at = 100 keV, as Fig. 1 in|Padovani et al.| (|2009[ ) shows. 
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15 1. The values used here are shown in Table 



tions 

[TJ The calculation of the ionization rate is discussed here at 
the example of W51C, but for the other objects the same 

procedure is done to obtain the result. 

In the case of W51C, for d = 6 kpc ( |Abdo et al.][2009| , 
E'lb = 1 GeV, a = 2.0 and E'^in = lO'^ eV, only the total 
proton energy budget of protons with a minimum kinetic 
energy of i?min, Wp and the average hydrogen density of the 
cloud are free parameters (see equation 15 1. The product 
of these two quantities needs to be 

Wpuu = 7.7 X 10^^ erg cm^^ (20) 

to produce the modeled gamma spectrum shown in Fig. [2] 
As can be seen, the modeled spectrum matches the detec- 
tions well. 

Because there is no precise estimate of the average den- 
sity of the cloud, here nn = 100 cm~'^ is assumed, thus 
offering the required value of Wp. The result for Wp sim- 
ply scales inversely with nu, if «h should turn out to be 
different. 

The normalization factor of the proton SED, Op is cal- 



culated following equation ( 14 ), where Wp is the total inter 



acting proton energy budget and the volume of the source 
K;ioud is assu med to be spherical with a radius taken from 
obser vations ( |Abdo et al.|2009l|2010c|a|d|b[ [Castro fc Slane 
20101. The lower integration limit indicates the minimum 
energy of particles contributing to ionization processes. As 
a conservative approximation, here i?min — 10 MeV is 
used. The result for each object is shown in Table [2] 

With this normalization one can compute the ioniza- 
tion rate du e to primary proton ionization performing the 
integration (17 1 for different values of E\\^ and a fixed value 
of ii'max — iGeV. The lower break energy is of large im- 
portance for the result, because the ionization cross section 
is the largest at low energies and declines rapidly with in- 
creasing energy. The upper integration limit may be any 
value from 1 GeV or higher because of the low ionization 
cross section for high energies. Changing -Emax to 1 PeV 
does not change the ionization rate significantly. However, 
changing Eib from 1 GeV to a different value results in a 
different ga mma spectrum normalization , as can be seen 
in equation ( |15[ ). 

Figure [3] sEows that even choosing a large value for the 
lower break energy would result in an ionization rate at 
least an order of magnitude greater than the Galactic aver- 
age form£kculBr_clouds^ reported to be about 2 x 10" 

by 



Neufeld et al. 



(2010), for at least two objects, namely 
W49B and 3C 391. TEe different lines refer to different val- 
ues for the power law index a. As can be seen, the value 
of a does influence the ionization rate, but only to a rather 
small extent. The ionization rate is much more sensitive to 
the choice of the lower break energy E\i,. If E\i, = 100 McV 
can be used, only the ionization rate for W51C would be 
lower than the Galactic average for molecular clouds, while 
the other objects would exceed this value for all spectral in- 
dices a. If even a value for as low as 30 MeV is suitable, 
the ionization rate for all objects would be greater than the 
Galactic average for molecular clouds. 

If even protons with a minimum energy of i?min = 
2 MeV could penetrate the cloud and thus contribute to 
the ionization, the ionization rates would increase further 



According to Indriolo et al. ( 2009 ) , this might be reasonable 
and would increase the ionization rate by 4 - 40% for a = 1 
and Eib = 30 MeV. For larger values of a and i?ib, this 
enhancement of the ionization rate would be significantly 
lower. 

The resulting ionization rates for each SNR known to 
be associated with a molecular cloud are given in Table [3j 

5. Uncertainties 

As mentioned above, since there is no observational data 
concerning the primary proton SED for energies below 
E ^ 1 GeV, the spectral behavior in this low energy 
domain is unknown. However, assuming a positive spectral 
index a = 2 is rather conservative and should offer a lower 
limit on the spectrum and therefore on the ionization rate. 
This is due to the fact that for an injection spectrum of 
oc p"**, for a loss term of p oc like ionization losses, 
would be modified as oc p^~'^, and the power law indices 
for the sources discussed are s = 1.5 - 2.45 . However, this 
is likely to be the largest uncertainty. Furthermore, there is 
no precise estimate of the average density of the molecular 
clouds, so a value of tih = 100 cm^'^ is assumed. This 
is also a major source of uncertainty. The primary proton 
SED scales linearly with the inverse value of rin (see equa- 
tion 12 1. Should the average density differ from the assumed 



value, then most likely it will be larger and therefore the pri- 
mary proton SED and the ionization rate would decrease. 
Another critical parameter is the volume of the cloud. This 
is a particularly difficult aspect, since only the primary pro- 
ton SED does depend on it (see equation 14 1, but not the 
resulting gamma spectrum (see equation 15 KThe radii used 
here refer to the whole SNRs and should therefore offer up- 
per limits on the cloud volume, which would result in lower 
limits for the primary proton SEDs and thus the ionization 
rates. With this respect, our calculations therefore repre- 
sent a conservative approach. The distance of the object 
from Earth enters the primary proton SED and the ioniza- 



tion rate as Ur, 



oc 



Earth— source 



(see equation 12 1, so the 



results are more sensitive on the lower break energy Ei^ 
than on the distance. 

Taking all the uncertainties discussed into considera- 
tion, still at least two sources, namely W49B and 3C 391, 
remain with an unusually large ionization rate at least one 
order of magnitude greater than the Galactic average for 
molecular clouds. Quite a few sources drop below Galactic 
average values in the most conservative scenario. This is 
not expected, since ionization at the accelerator itself is 
not likely to be lower than on average, but rather higher. 
Thus, we would rather expect the primary spectrum to ex- 
tend toward lower energies or the interaction volume to be 
significantly smaller. 



6. Signatures 

The enhanced ionization rate in the molecular clouds 
triggers a chemical network, forming a variety of ion- 
ized molecules in ro tationally and vibration ally excited 
states, see e.g. Black (1998 ) and Black (20071. Once these 
molecules are sufficiently abundant, their relaxation results 
in characteristic line emission. Though in principle similar 
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Fig. 2. Modeled gamma ray spectrum and Fermi-LAT observational data. The modeled spectrum shown is for 
E'lb = 1 GcV and a = 2.0, but the spectra for i?ib down to 30 MeV and a down to 1.0 practically coincide 
with the spectrum shown due to the low cross section below 1 GeV. 
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Table 2. Proton SED normalization Up for each source for different lower break energies i?ib and spectral indices a below 
the lower break energy in [erg~^ s~^ cm~^]. 
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Fig. 3. Ionization rates versus lower break energy Sib for different spectral indices below this break: a = 2.0 (solid 
black line), a = 1.0 (dotted red line). 
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Table 3. C^^/Cgai aver objccts, calculated for different spectral indices a below the lower spectral break E'lb and 

different lower break energies ii^ib. 



chemical signatures could be produced by X-ray ionization, 



X-rays have a much shorter penetration depth ( Beskin & 



Division 2003) and therefore are not capable of ionizing 
the cloud's interior as effectively as cosmic rays. This im- 
plies that the detection of chemical signatures similar to 



those shown in Becker et al. (201X1) would be a strong ar- 
gument for a proton SED capable of producing these large 
ionization rates, which in turn would be a very strong argu- 
ment for hadronic interactions as the dominant process in 
forming the detected gamma radiation in the GeV regime. 
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Uncertainties do not allow the prediction of the exact value, 
but the statement that in general, an enhanced ioniza- 
tion rate is expected. Observations of the objects with e.g. 
Herschel and ALMA will give a better idea of exact values 
in the future. First results of such observations are summa- 
rized in section [T] 

In steady state the number densities of the transient 
ions 0+, 0H+, H^O"*", and HeH"*" are exp ected to be pro- 
portional to ( Herbst fc Klemperer|l973 ). In fully molec- 
ular regions, where 



n(H) « n(H2) , 
of a low ionization level 

a;(e-) < 10"^ , 



(21) 



(22) 



competing processes almost do not interrupt the sequence 
of H-atom abstraction reactions. So nearly each ionization 
of hydrogen leads to the formation of Hj, H'^ , 011+ , H2O+, 
and HsO"*". The time-scales to achieve steady state are very 
short, 

-3 



1 - 1000 



cm 



n(H2) 



yr 



(23) 



These time-scales are shorter than the age of the SNRs 
dealt with here. 



As stated also in Becker et al. (2011 1, the concrete ion- 
ization signatures to expect from the enhanced ionization 
rates cannot be predicted quantitatively yet. However, the 
statement that enhanced ionization signatures are to be 
expected can clearly be made. These signatures would be 
characteristic rotation-vibration emission lines from abun- 



dant i onized molecules as e.g. and (see [Becker et al. 
(2011) for an example of such a spectrum). I'he environ- 



ments of SNR-MC systems are therefore considered to be 
optimal for a first time direct detection of H^. Should such 
ionization signatures be detected in spatial correlation with 
the GeV gamma ray emission from an SNR associated with 
a molecular cloud, this would provide a strong hint at in- 
elastic proton-proton scattering as the dominant process in 
forming these gamma rays. 



7. First experimental evidence 



In two sight lines in the IC443 complex, Indriolo et al 
(2010) found a large column density of Hj^, A^(H|^) = 
3 X lO^'* cm~^, indicating an enhanced ionization rate of 
^ 2 X 10^^^ s^^. This favors a less conservative cal- 
culation of the ionization rate in this region, as mentioned 
in section [5j Dense molecular gas associated with the W28 
SNR shows evidence of heating and chemistry driven by 
shock waves ( [Nicholas et al.| 2011[ [2012), but the high exci- 
tation of ammonia molecules m one molecular core can also 
be interpreted in terms of ion chemistry driven by an en- 
hanced rate of ionization as outlined in the next subsection. 
A derivation of the ionization ra te is given in subsection|7.1| 
In a molecular cloud of W51C, Ceccarelli et al. (2011) de- 
rived an enhanced ionization rate by measurements of the 
DCO+/HCO+ ratio of C"' ~ 10"^^ s'^, which also fa- 
vors a less conservative calculation of the ionization rate. 
These detections encourage additional observations toward 
the direction of the discussed SNR-MC systems, in partic- 
ular toward the very promising candidates W49B and 3C 
391. 



7.1. RADEX modeling for W28 

The W28 SNR is a prominent example of an association 
of partly resolved HESS 7-ray sources and dense, shocked 
molecular g revealed by the molecular line observa- 

tions of Nicholas et al. ( 2011[ 2012). In particular, the cm- 
wave inversion transitions of ammonia (NH3) show peak in- 
tensities that coincide with the positions of peaks in 7-ray 
emission and radio synchrotron emission. Conventionally 
NH3 emission is thought to probe dense molecular gas 
and the relative intensities of the inversion transitions are 
considered to be goo d indicators of kinetic temperature. 
Nicholas et al. ( |2011 ) identified a molecular cloud Core 2 
that is associated with the 7-ray source HESS J1801 — 233 
and where the {J,K) = (3,3) inversion line of NH3 is un- 
usually intense compared with the (1,1) and (2,2) lines 
of lower excitation. Moreover, they clearly detected emis- 
sion in the (6, 6) line and weakly in the (9, 9) line at the 
same position. From a three-dimensional radiative trans- 
fer an a lysis of their most sensitive NH3 spectra, |Nicholas| 
et al. (2011) determined a best-fitting hydrogen number 



density of nn = 10^ ''^ cm~'^ and a kinetic temperature 
T = 95 K, but excluded the weak (9, 9) line from the anal- 
ysis and noted that the model underestimates the intensity 
of the (6, 6) line. This analysis makes the standard assump- 
tion that the excitation is controlled by inelastic collisions 
of H2 with NH3 in competition with the radiative transi- 
tions. The NH3 lines in Core 2 have unusually large line 
widths, indicating Doppler velocity dispersions exceeding 5 
km s~^. Both the relatively high kinetic temperature and 
large velocity dispersion might result from the dynamical 
interaction of the expanding SNR with a dense molecular 
cloud. We have considered an alternative explanation of the 
high excitation of NH3, which might apply in a region of 
enhanced cosmic-ray ionization rate. A chemical formation 
source for the highly excited inversion levels would natu- 
rally account for the superthermal line widths: the newly 
formed molecules would be translationally hot - they gain 
kinetic energy from the enthalpy change in the chemical 
formation process. If NH3 is formed in a sequence of exoer- 
gic ion-molecule reactions that are initiated by cosmic-ray 
ionization of hydrogen and/or nitrogen, then the formation 
process itself, mainly 



nh:^ 



NH3 + H 



leaves the product NH3 molecules initially in highly excited 
rotational states, which relax rapidly by submm-wave rota- 
tional transitions. However, all steps in this radiative cas- 
cade that pass through rotational states (J, K) with J = K, 
leave molecules stranded in the lower inversion level, be- 
cause all these levels are highly metastable. If the rate of 
formation of NH3 is high enough, compared with the rates 
of inelastic collisions, then the formation process itself can 
account for observable populations of highly excited states. 
Indeed, this mechanism of "formation pumping" can mimic 
a coUisionally excited component of molecular gas at a tem- 
perature > 100 K, which would otherwise be needed to 
populate such metastable states as (6,6) and (9,9), which 
lie at energies E{J,K)/k = 408 K and 853 K above the 
ground state, respectively. 

We have calculated models of the non-LTE excitation 
of NH3 in which the formation and destruction processes 
are included together with the inelastic collisions involv- 
ing H2 and and all relevant radiative processes, through 
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use of the RADEX code as described by [van der Tak et aL] 
(2007). There are too many free parameters and too many 



uncertainties in colhsional rates to permit a fully optimized 
model. However, the observed intensities of the inversion 
lines of NH3 are reproduced fairly well in a model at ki- 
netic temperature Tk = 80 K with densities of H2, e^, 
and of 1000 cm^"^, 0.1 cm~^, and 0.03 cm^"^, respec- 
tively. The fractional abundance of ammonia relative to H2 
is 7.7x 10~* over a region of path length L = 1.3x 10^^ cm = 
4.2 pc, which corresponds to the extent of the strong emis- 
sion observed in the (3,3) inversion line in cloud Core 2. 
The Hj]' ions are included in order to account for the rel- 
ative amounts of NH3 in ortho-symmetry states, (3,3), 
(6, 6), (9, 9), and in para-symmetry states, (1, 1) and (2, 2), 
through reactive collisions that change the nuclear-spin 
symmetry of the molecule. In this model, the highly excited 
states are populated largely by the formation process itself, 
and the inferred rate of formation is 7.7 x 10"^'^ ammonia 
molecules cm"'^ s~^. This is balanced by a destruction rate 
of approximately 10~* s~^. These rates would require that 
C^^/C^l ~ 100 and thus imply a lower spectral break in 
the cosmic ray spectrum at En, 10 MeV (see Table |3]). 

The high abundance and high excitation of NH 3 ob- 
served in molecular Core 2 by Nicholas et al. ( |2011 ) could 
be explained as the result of rapid formation in an ion- 
driven chemistry. As discussed by Nicholas et al., the NH3 
might result from shock-driven chemistry, which would also 
fit well with the other tracers of molecular shocks that they 
describe. The clearest way to distinguish between these two 
explanations would be observations of molecular ions like 
Hj^ and H3O+, which would be expected to have greatly 
enhanced abundances if the cosmic ray ionization rate is 
unusually high. 



8. Conclusions and outlook 

In this paper we compute ionization rates for molecular 
clouds known to be associated with SNRs based on the as- 
sumption that the GeV gamma ray emission from these ob- 
jects is due to neutral pion decay formed in proton-proton 
scattering in the molecular clouds. The computed ioniza- 
tion rates for at least two objects are above Galactic aver- 
age for molecular clouds in the most conservative scenario. 
Therefore ionization signatures in the form of rotation- 
vibration line emission from molecular ions are likely to 
be detected from these two objects. The spatial correla- 
tion of the detection of GeV gamma rays on the one hand 
and rotation- vibration line emission from molecular ions on 
the other hand would strongly hint at the hadronic origin of 
the detected GeV gamma ray excess, explaining both detec- 
tions by one population of cosmic ray particles. However, 
there is the caveat that low energy CRs are efficient in 
ionizing, whereas high energy CRs are responsible for the 
gamma radiation. One can expect that low and high energy 
CRs are accelerated in the same objects, but an extrapo- 
lation of the high energy CR spectrum inferred from the 
gamma ray detections to low energies is not exempt from 
problems. Still, recent observations hint at enhanced ion- 
ization rates in molecular c louds associated with S NRs, e.g. 
Nicholas et al. (2011) and Ceccarelli et al. (2011), in sup- 
port of the presented model 

In future work, differential propagation of the primary 
protons into the molecular cloud as well as the considera- 



tion of secondary ionization are planned to be implemented 
to offer more precise predictions concerning the ionization 
signatures to be expected. 
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